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We examine a simple model of Pb(Ini/2Nbi/2)03 (PIN), which includes both long-range 
dipole-dipole interaction and random local anisotropy. An improved algorithm optimized for 
long-range interaction has been applied to efficient large-scale Monte Carlo simulation. We 
demonstrate that the phase diagram of PIN is qualitatively reproduced by this minimum model. 
Some characteristic features of relaxors such as nanoscale domain formation, slow dynamics, 
and dispersive dielectric responses are also examined. 
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Relaxors, whose first discovery dates back to over half 
a century ago,^^ have newly attracted much interest over 
the past few decades owing to their colossal dielectric 
and piezoelectric responses that are appealing for indus- 
trial applications. Despite intensive research, however, 
the physical origin of the unusual properties of relaxors 
is not fully understood yet.^^ This is because nanoscale 
intrinsic randomness in relaxor systems has to be ad- 
dressed appropriately. 

As a simple example to understand the physics 
of relaxors, let us choose the lead-based relaxor 
Pb(Ini/2Nbi/2)03 (PIN), which is the target of the 
present theoretical work. The intrinsic inhomogeneities 
of relaxor systems originate from the random configura- 
tion of In^+ and Nb^+ at the B-site in the perovskite 
structure, while the dielectric property of such systems 
is governed mainly by the replacement of Pb and O ions. 
This system has the advantage that the strength of ran- 
domness can be controlled by adjusting annealing tem- 
perature.^^ For the lowest annealing temperature, an al- 
ternate order of In and Nb atoms stabilizes the fourfold 
antiferroelectric (AFE) phase. For the higher annealing 
temperature, an increase in B-site randomness decreases 
AFE transition temperature. Under sufficiently strong 
randomness in the B-site, the ferroelectric phase accom- 
panies relaxor properties. This behavior of PIN suggests 
that the potential FE instability exists behind the AFE 
phase, and that its emergence by the selective suppres- 
sion of the AFE phase due to B-site randomness is rele- 
vant to FE relaxors. Recent experimental progress in the 
X-ray diffraction analysis facilitated by strong-intensity 
photon sources has also enabled us to detect the poten- 
tial competition between the FE and AFE phases of PIN 
from the viewpoint of lattice dynamics.^' 

There are several theoretical methods for elucidat- 
ing the emergence of relaxor properties. They uti- 
lize the existing effective theory dealing with a strong 
randomness such as the extension of the Ginzburg- 
Landau-Devonshire theory,^' the spherical random- 
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bond-random-field model, and the dynamical-related 
model. ^'^^^ Although these approaches have provided 
useful description of relaxor properties, it is unclear 
how their model parameters should be derived micro- 
scopically. At present, first-principles calculation is play- 
ing a crucial role in the microscopic description of fer- 
roelectrics and ferroelectric relaxors. In the pioneer- 
ing work by Cohen and coworker, the electron- 
structure calculation of BaTiOa and PbTiOa indicated 
the importance of the covalent nature among composite 
atoms through Ti-0 hybridization. After the success of 
their work, microscopic first-principles approaches have 
been adopted in a number of issues. Further theoretical 
progress has been achieved by a hybrid method com- 
bining the microscopic derivation of the effective Hamil- 
tonian and its numerical simulation. Zhong et al. have 
derived the effective Hamiltonian of BaTiOs by extract- 
ing the adiabatic potential of important phonon modes 
from first-principles calculation.^^' Monte Carlo simu- 
lation for this effective Hamiltonian has enabled the suc- 
cessful reproduction of the sequential finite-temperature 
phase transitions of BaTiOa in good agreement with ex- 
perimental results. Recently, the hybrid method has also 
been applied to the study of relaxor systems. 

Thus, the utilization of first-principles calculation is a 
promising method for understanding relaxors. However, 
there still remain several difficulties to overcome for the 
actual application of this method to the study of relaxors. 
One major difficulty is in the use of the numerical solver 
of the effective Hamiltonian for phonons. Even though 
the effective Hamiltonian derived for relaxors is in a sim- 
ple form, its numerical simulation may become extremely 
difficult because it is inevitable to encounter slow dynam- 
ics inherent to random systems. The Monte Carlo sim- 
ulation of the effective Hamiltonian needs a large num- 
ber of Monte Carlo steps owing to its long correlation 
time. Therefore, the development of appropriate numeri- 
cal solvers for relaxor systems seems indispensable. This 
situation may remind us of Monte Carlo studies of spin- 
glass systems. Although the Hamiltonian of spin-glass 
systems has a simple form, much effort has to be exerted 



2 J. Phys. Soc. Jpn. 



Letter 



Yusuke Tomita et al. 



for sufficient Monte Carlo sampling to understand their 
peculiar slow dynamics. In the research field of relaxors, 
however, the importance of developing a numerical solver 
has not been addressed so far. 

In this paper, we propose a simple model of ferroelec- 
tric relaxors to explain the phase diagram of PIN. Our 
model includes both long-range dipole-dipole interaction 
and local randomness. We apply the new effective algo- 
rithm optimized for long-range interaction proposed by 
Fukui and Todo.^^^ This new algorithm enables us to 
simulate long-range interaction systems with the cost of 
O(A^) with respect to the system size TV, while the con- 
ventional simulation for the same system takes the cost of 
0(A^^). In addition to the substantial reduction in com- 
putational time by employing this algorithm, we adopt 
the exchange Monte Carlo method to attenuate the slow 
dynamics of random systems. We show that our sim- 
ple model may qualitatively reproduce the phase diagram 
of PIN. Moreover, we demonstrate that some character- 
istic features of relaxors such as domain formation and 
dispersive dielectric response are reproduced reasonably. 
Our results indicate the power of the sophisticated Monte 
Carlo method as well as the possibility that relaxor sys- 
tems may be represented by a simple model. 

We consider the model Hamiltonian for PIN on a 2D 
square lattice as 
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Fig. 1. (Color online) (a) A plot of four-fold antiferroelec- 
tric(AFE) order parameters in the absence of randomness {p = 0) 
and (b) a plot of ferroelectric(FE) order parameters at p = 0.5 
are shown for different system sizes N = L x L. (c) Phase dia- 
gram obtained by the minimum model. Phase boundaries in the 
figure are visual guides. 



where Si is a unit vector in the x?/-plane representing 
the dipole moment on the ith unit cell induced by the 
off-center replacement of the Pb atom. The first term 
of the Hamiltonian is the dipole-dipole interaction de- 
pendent on the relative position Vij = rj — Vi between 
the sites i and j, while the second term describes local 
anisotropy whose direction and strength are denoted as 
Di. In order to reproduce the phase diagram of PIN, we 
design our model such that the FE phase is stabilized 
by dipole-dipole interaction, while the AFE phase is sta- 
bilized in the presence of an alternative change in Di. 
Supposing that B-site randomness affects only the local 
energy change through the anisotropy parameter Di^ we 
can expect that all the features of PIN, i.e., the AFE 
transition in the ordered PIN, its suppression by B-site 
randomness, and the appearance of the FE domain in 
the sufficiently disordered PIN are reproduced. 

The detailed setting of our model is given as follows. 
In a naive 2D square lattice, the dipole-dipole interac- 
tion does not lead to ferroelectric instability because 
the columnar antiferroelectric state is more favored in 
orthogonal lattices. Therefore, we modify our model 
slightly: We divide the unit cells into two bipartite sub- 
lattices named P and Q, and shift only the Q-sites in the 
z-direction by a unit length. This rearrangement of unit 
cells ensures ferroelectric instability in the absence of lo- 
cal anisotropy. The antiferroelectric instability is driven 
by an alternative arrangement of two types of anisotropy. 



namely, Di and D2 in the P- and Q-sites, where Di = 
{D/V2,-D/V2) and D2 = {D/V2,D/V2). The ran- 
domness of D is controlled by the probability p, where 
Di {D2) is attached to the P-(Q-)site with a probability 
(1 — p)^, oppositely to p^, and is turned off (D=0) with 
2p{l —p). The two limiting values p = and p = 1/2 cor- 
respond to the ordered and completely disordered PINs, 
respectively. 

The model Hamiltonian is examined by Monte Carlo 
simulations under a periodic boundary condition to min- 
imize the strong surface effect due to long-range inter- 
action. Long-range interactions are summed up by the 
Ewald summation technique. We employ the 0{N) 
method based on Walker's algorithm for efficient up- 
date^^'^^^ as well as the temperature-exchange algo- 
rithm. One exchange trial between replicas was made 
for each of the 10 MC steps. The system size is taken up 
to AT = 32 X 32, for which 3 x 10^ MC steps for ther- 
malization and 2 x 10^ MC steps for measurement were 
needed in the severest case of p = 1/2. The sample av- 
erage is taken over 10 different random configurations of 
Di for each p. Throughout this paper, the strength of 
the anisotropy is fixed as \Di\ = 1. 

Figure 1(a) shows the temperature dependence of the 
average squared four-fold staggered polarization ttiaf for 
the ordered PIN [p = 0). The pattern of the AFE order- 
ing, which is shown in Fig. 1(c), agrees with experimental 
observations. As the system size N = L x L increases. 
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Fig. 2. (Color online) Log-log plot of the normalized time auto- 
correlation functions. While the autocorrelation functions decay 
smoothly for the no-anisotropy system at T = 0.2 and the ran- 
dom system at T = 0.3, the autocorrelation function shows a 
slow decay for the random system at T = 0.2. 



a sharp increase in the order parameter appears below 
a critical temperature, indicating the phase transition. 
The transition point is determined by the crossing point 
of the Binder parameters, (^af)/(^af)^' of L = 16 
and 32. By a similar analysis of nonzero values of p, 
the phase boundary of AFE in the p-T plane is deter- 
mined, as shown by full circles in Fig. 1(c). For small 
values of p, the transition temperature of AFE is sup- 
pressed by the B-site randomness. For sufficiently large 
values of p, the FE domain develops at low temperatures 
instead of the AFE phase. Figure 1(b) shows the average 
squared uniform polarization as a function of temper- 
ature in the completely disordered case {p — 1/2). An 
abrupt increase in squared polarization below a thresh- 
old temperature for the = 32 x 32 system indicates a 
rapid development of the FE domain. In our simulation, 
however, no long-range FE ordering could be detected 
by finite-size scaling because of the very slow relaxation 
of the Monte Carlo sampling and the complex size de- 
pendence in the presence of strong randomness. Here, we 
roughly estimate the threshold temperature at which the 
local FE instability rises up by Binder-parameter analy- 
sis between L = 16 and 32, and plot it using the empty 
circle in Fig. 1(c). At approximately p = 0.4, the com- 
petition between AFE and FE makes the Monte Carlo 
calculation severe, and the phase boundary could not 
be determined. Rough interpolations for AFE and FE 
are shown by the solid curves in Fig. 1(c) only as visual 
guides. The obtained phase diagram is in good agreement 
with experimental results of PIN. 

For large values of p, some characteristic features of 
relaxor systems appear below a crossover temperature, 
which locates higher than the threshold temperature 
for FE domain formation, as indicated roughly by the 
dashed line in Fig. 1(c). The emergence of relaxor prop- 
erties is indicated significantly by the slow convergence of 
the Monte Carlo average. Figure 2 shows the normalized 
autocorrelation of the energy £^ for A/" = 32 x 32 defined 
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Fig. 3. (Color online) Snapshots of dipole configuration for p = 
0.5 and T = 0.1. Each color denotes the angle of dipole moment. 
Reference color bars are shown on the right. The bottom color 
denotes — tt, and the top color denotes tt. The last graph shows 
the result for the same parameters but for a model with a dipole 
interaction interrupted up to the next-nearest sites. 
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as a function of the Monte Carlo step t in the completely 
random case {p = 1/2). Here, for the purpose of exam- 
ining the dynamics of the system, we turned off the ex- 
change process between replicas. We find that the de- 
cay of the autocorrelation function is extremely slow at 
T = 0.2 in clear contrast with that at a higher temper- 
ature T = 0.3. This slow energy relaxation indicates a 
glasslike behavior for relaxor s at low temperatures. In 
order to confirm that this slow relaxation comes from 
randomness of the B-site, we also calculated the auto- 
correlation at T = 0.2 in the absence of local anisotropy 
{\Di\ = 0), for which normal FE ordering is expected. As 
seen in Fig. 2, the autocorrelation function follows a sim- 
ple decay with no anomalous slow relaxation. Thus, the 
present model is expected to exhibit relaxor properties 
under sufficiently strong randomness below a crossover 
temperature. 

In order to visualize the glassy state realized in the 
relaxor phase, we show in Fig. 3 the spatial pattern of 
dipole directions by taking snapshots as a function of 
the MC steps for p = 0.5 and N = 32 x 32. We can 
see that mosaic-like FE domains are formed after a finite 
relaxation time. The boundaries (domain-wall regions) 
between the neighboring FE domains are rather clear. 
One may consider that the FE domain structure is deter- 
mined by a partially ordered configuration of anisotropy, 
the so-called chemical nanoregion (CNR). The present 
FE domain is, however, irrelevant to CNR, because par- 
tial ordering is absent for a complete random choice of 
anisotropy at each site in our model. 

As has already been mentioned, the present system 
possesses FE instability in the absence of local anisotropy 
{D = 0). The random configuration of anisotropy simi- 
larly causes a local FE instability due to effective cancel- 
lation of anisotropy, but at the same time prevents the 
growth of an FE domain into a uniform FE order, as 
seen in Fig. 3. A small but finite change of the snapshot 
in the long-time region of the Monte Carlo simulation 
indicates that a very slow fiuctuation of the FE domain 
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governs the relaxor properties in the present model. Here, 
we should note that the evolution of the dipole configu- 
ration along the Monte Carlo step is not directly related 
to actual real-time dynamics. The snapshots, however, 
provide intuitive and fruitful insights for understanding 
physical processes in relaxor systems. In our model, the 
long-range nature of the dipole interaction is essential 
for the formation of the FE domain structure. To see 
this, we show a snapshot of Monte Carlo simulation for 
a model with dipole interaction interrupted up to the 
next-nearest site in the last graph of Fig. 3. We find that 
the domain structure disappears in this modified model. 
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Fig. 4. (Color online) Dielectric constants evaluated as functions 
of temperature under the ac electric field for several frequencies 
in the (a) ordered case (p = 0) and (b) completely disordered 
case {p = 1/2). 



The existence of nanoscale frustrated FE domains in- 
dicates the high degeneracy of glassy states. This degen- 
eracy is expected to be responsible for large dielectric 
responses since a significant change in the polarization 
of such FE domains is possible by a small external elec- 
tric field. To see this, we measure the Monte Carlo evo- 
lution of the total polarization under a small 'ac' electric 
field, which varies periodically along the Monte Carlo 
step. We then calculate the 'ac' dielectric constant in 
the linear response regime. We should note again that 
the 'ac' dielectric constant thus obtained is not equal 
to the real-time ac dielectric constant. Furthermore, we 
should keep in mind that this is a nonequilibrium re- 



sponse since the relaxation to equilibrium states is never 
realized in the glassy phase in the present simulation. 
Nevertheless, the 'ac' dielectric response calculated here 
is expected to mimic the actual ac response in a qualita- 
tive level. In Fig. 4, we show the 'ac' dielectric constant 
as a function of temperature for three frequencies in the 
ordered {p — 0.0) and completely disordered {p = 1/2) 
cases. In both cases, dielectric constant shows a peak 
at approximately the transition (crossover) temperature. 
The marked difference is found in the low-temperature 
phase. In the ordered case, dielectric constant sharply 
drops in the low-temperature AFE phase, and becomes 
almost independent of frequency. On the other hand, in 
the disordered case, it decreases gradually as temper- 
ature decreases, and a strong frequency dependence re- 
mains. This strong dispersion of dielectric constant at low 
temperatures can be explained as follows. Under a high- 
frequency electric field, each dipole inside FE domains 
may respond to an external field. Under a low-frequency 
electric field, however, those with frustration that con- 
struct large domains start to respond. Therefore, in a 
minimum model, a dipolar glass with ferroelectric order- 
ing is realized. This ordering causes a broadened dielec- 
tric constant and a strong dependence on frequency in 
the relaxor phase. 

In summary, we examined a simple theoretical model 
of PIN composed of dipolar interaction and local ran- 
dom anisotropy. We demonstrated that efficient Monte 
Carlo simulations equipped with an improved algorithm 
optimized for long-range interaction may access several 
characteristic features of relaxor s. The phase diagram of 
PIN was qualitatively reproduced by appropriate inclu- 
sion of the intrinsic competition between the AFE and 
FE phases. By the examination of the Monte Carlo evo- 
lution of the dipole configuration, we demonstrated some 
of the glassy behaviors inherent to relaxor systems such 
as the FE domain formation, extremely slow dynamics, 
and strong frequency dependence of dielectric responses. 

We end this paper by mentioning the future outlook of 
theoretical approaches to studying relaxors. The model 
we treated here includes a few artificial assumptions. 
They, however, will be removed by replacing our model 
with the effective Hamiltonian derived by first-principles 
calculation. We stress that the smart Monte Carlo al- 
gorithm is applicable not only to the rotator model but 
also to the continuous- variable model. The hybridization 
of the first-principles calculation and statistical approach 
based on the Monte Carlo simulation will be an effective 
means of elucidating the microscopic origin of relaxors. 
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